eststo clear


if("`c(os)'"=="Windows"){
	local tableSaveDir = "ReviewResponseResults\"
	local plotSaveDir = "ReviewResponseResults\"
	local pythonSaveDir = "PythonScripts\"
}
else{
	local tableSaveDir = "ReviewResponseResults/"
	local plotSaveDir = "ReviewResponseResults/"
	local pythonSaveDir = "PythonScripts/"
}

//Establish Control Variable Sets
//Spec1 controls are the same for both registrations and VMT 
local spec1Controls = "c.normYear#i.stateGroup c.quadYear#i.stateGroup asinhMetroPop asinhNonMetroPop asinhMetroInc asinhNonMetroInc logrealmeangasprice logemployment loglicenseddrivers logrealstategdp logroadmileage"

//Other controls are not the same for both registrations and VMT 
//Registration Controls
local spec2RegControls = "`spec1Controls' i.transactionDataDummy"
local spec3RegControls = "`spec2RegControls' i.regDummyExtSpec3"
local spec4RegControls = "`spec2RegControls' i.regDummyExtSpec4"


di "`spec3RegControls'"

//VMT Controls
local spec2VMTControls = "`spec1Controls'"
local spec3VMTControls = "`spec2VMTControls' i.vmtDummyExtSpec3"
local spec4VMTControls = "`spec2VMTControls' i.vmtDummyExtSpec4"

//Absorb Vars are always time and state fixed effects
local absorbVars = "ib51.stateGroup ib1995.year"

//Establish Cluster Variable
local clusterVar = "stateGroup year"

//Establish primary null
local primaryNull = 1.6
local altNull = 1
local altNull2 = 0.3
local altNull3 = 0
local altNull4 = -2.5

frame create ruralUrban
frame change ruralUrban

if("`c(os)'"=="Windows"){
	local dataFile = "CleanedData\UrbanRuralVMT.csv"
}
else{
	local dataFile = "CleanedData/UrbanRuralVMT.csv"
}

import delimited using "`dataFile'"

frame change default
frlink 1:1 year state, frame(ruralUrban)
frget ruralvmt urbanvmt, from(ruralUrban)
gen logruralvmt = log(ruralvmt)
gen logurbanvmt = log(urbanvmt)


//Run urban model to get appropriate N and K values.
qui reghdfe logurbanvmt nosafetyind `spec2VMTControls', absorb(`absorbVars') cluster(`clusterVar')
local k_u = e(df_a)+e(df_m)
local N_u = e(N)
local r2_u = e(r2)
local N_clust_u = e(N_clust)
local correction_u = (`N_clust_u')/(`N_clust_u'-1)*(`N_u'-1)/(`N_u'-`k_u')

//Run rural model to get appropriate N and K values.
qui reghdfe logruralvmt nosafetyind `spec2VMTControls' if state!="Dist. of Col.", absorb(`absorbVars') cluster(`clusterVar')
local k_r = e(df_a)+e(df_m)
local N_r = e(N)
local r2_r = e(r2)
local N_clust_r = e(N_clust)
local correction_r = (`N_clust_r')/(`N_clust_r'-1)*(`N_r'-1)/(`N_r'-`k_r')

//Partial urban model variables
qui reghdfe logurbanvmt `spec2VMTControls', absorb(`absorbVars') cluster(`clusterVar') res(logurbanvmt_res)
qui reghdfe nosafetyind `spec2VMTControls', absorb(`absorbVars') cluster(`clusterVar') res(nosafetyind_res_u)

//Partial urban model variables
qui reghdfe logruralvmt `spec2VMTControls' if state!="Dist. of Col.", absorb(`absorbVars') cluster(`clusterVar') res(logruralvmt_res)
qui reghdfe nosafetyind `spec2VMTControls' if state!="Dist. of Col.", absorb(`absorbVars') cluster(`clusterVar') res(nosafetyind_res_r)

//Estimate partialled urban model
regress logurbanvmt_res nosafetyind_res_u, nocon
estimates store urban

//Estimate partialled rural model 
regress logruralvmt_res nosafetyind_res_r if state!="Dist. of Col.", nocon
estimates store rural

suest urban rural, vce(cluster stateGroup)
mat var1 = e(V)
suest urban rural, vce(cluster year)
mat var2 = e(V)
suest urban rural, vce(robust)
mat mat_beta = e(b)
mat varR = e(V)

mat mat_var = var1+var2-varR

matlist mat_beta
matlist mat_var

mat diff_r = [1, 0, -1, 0]
matlist diff_r

mat diff_value = diff_r*(mat_beta')
matlist diff_value

mat diff_var = diff_r*mat_var*(diff_r')
matlist diff_var

local chi_stat = (diff_value[1,1]^2)/sqrt(diff_var[1,1])

di "`chi_stat'"

local p_value = 1-chi2(1, `chi_stat')

di "`p_value'"
